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I. Introduction 


Ultracold atomic quantum gases in optical lattices have 
reached an unprecedented degree of control providing di¬ 
rect experimental access to a plethora of non-equilibrium 
phenomena n ■ This control includes the modula¬ 
tion of the interparticle interactions via confinement- 
induced, magnetic and optical Feshbach resonances [f§ 
the design of arbitrarily shaped optical traps with 
variable lattice depths, and the ability to move time- 
periodically or even accelerate the entire lattice struc¬ 
ture. This level of control and accuracy over the sys¬ 
tem parameters has opened the possibility to simulate 
and study quantum many-body phenomena in part in¬ 
spired from condensed matter physics 0 - For instance, 
when accelerating an optical lattice, representative pro¬ 
cesses are Bloch-oscillations Wannier-Stark lad¬ 

ders [TaGi, Landau-Zener tunneling 0,001 and photon 
assisted tunneling [0 ], to name only a few. A promising 
technique is the lattice shaking which has been used in or¬ 
der to address e.g., the coherent control of the superfluid 
to Mott insulator phase transition j2C| , parametric am¬ 
plification of matter waves 0|. four-wave mixing 0,0], 
topological states of matter (24| . hybridized band struc¬ 


ture 


25j, and even the engineering of artifici al g auge 


fields [26[. More recently it has been shown [27], |28| that 


one can use lattice shaking to probe coherent band cou¬ 
pling and realize the formation of ferromagnetic domains. 
Moreover, the dynamics induced by shaking an optical 
lattice can lead to an admixture of excited orbitals 0] 
and constitutes an emergent branch of modern quantum 
physics. 


A substantial part of the previous studies has been 
primarily focussing on the renormalization of the physics 
due to driving, the mean-field approach 0 for weak in¬ 
teractions, where the Gross-Pitaevskii equation is still 
valid, and a linear response treatment [30 ]. However, 
a relatively large modulation of the strength or of the 
frequency of the driving as well as strong interactions, 
calls for alternative methods which can take into account 
higher-orbitals. Indeed, the inclusion of higher-band con¬ 
tributions introduces new degrees of freedom and as a re¬ 
sult additional physical processes come into play. Hereby, 
a sinusoidal shaking of the optical lattice is a natural 
starting point which induces an in-phase dipole mode on 
each site. An interesting and so far largely unexplored 
direction is the study of the interplay between higher 
bands for the intra-well mode and the inter-well tunnel¬ 
ing dynamics with respect to the driving frequency, and 
the investigation of the effect of the interatomic interac¬ 
tions in the overall process. In this way, it is natural to 
start with the investigation of the few body analogue in 
order to achieve a more comprehensive understanding of 
the microscopic properties of the strongly driven inter¬ 
acting system. Although the major part of the presented 
results is devoted to the case of four bosons in a triple¬ 
well setup, we provide strong evidence that our findings 
are still applicable for larger lattice systems and larger 
particle numbers. 

Motivated by the recent experimental progress 0,0 
we investigate in the present work the effects a periodi¬ 
cally driven one-dimensional optical lattice can introduce 
in a small ensemble of ultracold bosons. The dynamical 
response of the system for a wide range of driving fre- 
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quencies is studied by means of the concept of fidelity or 
autocorrelation function. Even though we consider a sce¬ 
nario with a deep lattice such that the tunneling modes 
have a minor influence on the overall dynamics, a quite 
rich excitation spectrum is found. We note that such 
intra-band excitations, which lead to a coupling between 
the two lowest energy bands, have been exploited in order 
to realise single- and two-qubit gates, where the quantum 
bit has been encoded in the localised Wannier functions 
of the two lowest energy bands of each lattice site (3l|. In 
order to analyze the intra-well dynamics we employ the 
one-body reduced density matrix. The Fourier spectrum 
of the local one-body density as well as of the on-site den¬ 
sity oscillations are employed in order to obtain insights 
into the excited intra-well modes. We find a resonant 
behaviour of the dipole mode indicating that the intra¬ 
well dynamics can be controlled by adjusting the driving 
frequency. Moreover, the magnification of the intra-well 
generated mode at resonance is also manifested in the 
population of additional lattice momenta. Our investiga¬ 
tion of the resonances is supported by a Floquet analysis 
for the effective single-particle degree of freedom. This 
allows us to further explore the on-site dynamics and the 
inter-well tunneling that occur due to the driving. In¬ 
cluding interatomic interactions for larger atom numbers 
we analyze similarities and differences with respect to the 
single-particle description. The above outlined findings 
are confirmed for different filling factors, lattice poten¬ 
tials, and boundary conditions. To solve the underlying 
many-body Schrodinger equation we apply the ab—initio 
MultiConfiguration Time-Dependent Hartree method for 
Bosons (MCTDHB) [32| [ 33 } which is especially designed 
to treat the driven out-of-equilibrium quantum dynamics 
of interacting bosons. 

This article is organized as follows. In Sec.II we in¬ 
troduce our setup and the multi-band expansion. Sec. 
Ill contains the driven quantum dynamics first from a 
single-particle perspective, by performing a Floquet anal¬ 
ysis, and second by inspecting the dynamics of a small 
bosonic ensemble including repulsive interactions. We 
summarize our findings and provide an outlook in Sec.V. 
Appendix A briefly outlines our computational method. 

II. Hamiltonian and multi-band expansion 

This section is devoted to a brief presentation of the 
theoretical framework of our study. In particular, we 
shall briefly discuss the driven optical lattice, the under¬ 
lying many-body Hamiltonian, and the concept of multi¬ 
band expansion. The latter will be a useful tool in order 
to understand the excitations involved in the dynamics. 

A. Modeling the periodically-driven potential 

The periodic driving of an optical lattice can be ac¬ 
complished in two different ways. Retroreflecting mirrors 


that are used to form the lattice can be moved periodi¬ 
cally in space or, alternatively, a frequency difference be¬ 
tween counterpropagating laser beams can be induced by 
means of acousto-optical modulators [27j which renders 
the lattice time-dependent. Here, we model the driven 
optical lattice with a sinusoidal function of the form 

V 8h (x,t) = V 0 sm 2 [k 0 (x - Asinwpt)]. (1) 

Such a potential has been implemented in the experiment 
of e.g. ref. [2lj| . It is characterized by the barrier depth 
Vq, a lattice wave-vector k 0 = y, where l denotes the 
distance between successive potential minima, the am¬ 
plitude A and the frequency lod = 2tt/Td of the driving 
field. In an experiment ko is the wave vector of the laser 
beams which form the optical lattice, while its depth Vq 
can be tuned by adjusting the lasers intensity. 

B. The Hamiltonian 

The Hamiltonian of N identical ultracold bosons of 
mass M confined in a driven one-dimensional m-well op¬ 
tical lattice reads 

N Ti 2 d 2 

H = ^ — X j)> ( 2 ) 

i =1 1 i<j 

where Vi nt (xj — Xj ) = gm^ixi — Xj ) denotes the short- 
range contact interaction potential between particles lo¬ 
cated at position Xi, i = 1,2,...,TV. In the ultracold 
regime the interaction is well described by s-wave scat¬ 
tering whose effective ID coupling strength Q is given 

by gi D = f^(l - IC( X^° ) • Here “-l = V^T 

is the transverse harmonic oscillator length with wj_ the 
frequency of the two-dimensional confinement, while ao 
denotes the free space 3D s-wave scattering length. In 
this way, the interaction strength can be tuned either via 
ao with the aid of Feshbach resonances 9, 10], or via the 
transversal confinement frequency w_l 1341 l35j . 

For the sake of simplicity and computational conve¬ 
nience, we rescale the Hamiltonian (2) in units of the re- 
n 2 k 2 

coil energy Er = Then, the corresponding length, 

time and frequency scales are given in units of fcjj' 1 , 
ujft 1 = hEft 1 and ujr respectively. In our simulations we 
have used a sufficiently large lattice depth with values 
ranging from Vo = 4.5 Er to 8.0 Er such that each well 
includes three localized single-particle Wannier states. In 
particular, due to the deep optical lattice and small driv¬ 
ing amplitudes (in comparison to the lattice constant) 
mainly used in our simulations highly energetic excita¬ 
tions above the barrier are excluded and as a consequence 
heating processes can be minimized. The confinement of 
the bosons in the m-well system is imposed by the use of 
hard-wall boundary conditions at positions x a = 
where the potential is maximum. In addition, we set 
also Ti = M = ko = 1 and the coupling strength becomes 
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g = while A represents the dimensionless driving 
amplitude. The rescaled shaken triple well is given by 
V s h(xi,t) = Vosin 2 (a:j — Asinivijt) with the hard wall 
boundaries located at x a = ±37r/2. 

C. The multi-band expansion 

The understanding of the spatial localization of states 
in lattice systems makes the use of multi-band Wannier 
number states crucial as it includes the information of ex¬ 
cited bands and allows to interpret both intraband and 
interband processes. In general, this representation is 
valid when the lattice potential is deep enough such that 
the Wannier states between different wells have a very 
small overlap for not too high energetic excitation. In the 
present case where the potential is periodically driven the 
above description can still be used as long as the driving 
amplitude is small enough in comparison to the lattice 
constant l, i.e. A <C l. In this way, each localized Wan¬ 
nier function can be still adapted and assigned to a cer¬ 
tain well and the respective band-mixing is fairly small. 
For large displacements one should use a time-dependent 
Wannier basis in order to ensure that the correspond¬ 
ing on-site Wannier states are well-adapted to each well 
during the driving. 

To introduce the formalism, let us consider a system 
consisting of N bosons, m -wells and k localized single 
particle bands [36;, :37]. Then, the expansion of the many- 
body bosonic wavefunction in terms of the number states 
of non-interacting bosons reads 

|*} = E (3) 

W,I 

where \Ni, N 2 , ■■■, N m ) 1 is the multiband Wannier num¬ 
ber state and the element iV,; denotes the number of 
bosons being localized in the z-th well satisfying the con¬ 
straint X^=i Ni = N. The summation is performed over 
the different configurations of the N bosons according 
to their energetical order denoted by the index I. In 
particular, the index I corresponds to a high dimen¬ 
sional quantity I = (/ 1 , / 2 ,...,/ m ) which contains m el¬ 
ements each of them being a /c-component vector. More 
precisely, the g-th element can be written as I q =(I^\ 
Iq 2 \..., Iq k ■*), where I q ^ refers to the number of bosons 
located at the g-th well and fc-th band, satisfying the 

constraint J2 q Li J2i =1 E = N. Within the above no¬ 
tation one can investigate, among others, the probabil¬ 
ity of No < N bosons to be in an excited band or to 
find a specific number state configuration. Indeed, sup¬ 
pose the case of Nq < N bosons excited in the i -th band 
while the rest N — No lie in lower bands. Then, it must 
hold l[^ = = ... = Im = 0 for every j > i, while 

if 5 + if 1 + ... + itt = N 0 and if } + ... + if l} = TV — N 0 
for every j\ < i. 

Let us consider an example of a system with four 
bosons (TV = 4) confined in a triple well (m = 3) which 


includes three bands (k = 3). Then, for instance, the 
state 11,2, l)j with I = ( I l ,I m Jr ), and I L =I R =(0,1,0), 
iM=(0,l,l) denotes a state for which in the left (right) 
well one boson occupies the first excited band, whereas in 
the middle well one boson is localized in the first excited 
and one in the second excited band. As a final attempt, 
here, we make a link between the ground state and its 
dominant spatial configuration in terms of the aforemen¬ 
tioned multiband expansion. To do that, let us choose 
again a system consisting of four bosons in a triple well 
as it will extensively be used in the following. It is known 
that, in general, the ground state configuration depends 
on the interaction strength, while for the present system, 
i.e. N = 4 and m = 3, the on-site interaction effects 
will always be prominent. For the non-interacting case 
(g = 0) the dominant spatial configuration of the system 
is 11,2, l)j, with I L = I R = (1,0,0) and I M = (2,0,0) 
due to the hard-wall boundaries which render the middle 
and outer sites non-equivalent. In the course of increas¬ 
ing interaction a tendency towards a uniform population 
of each site, e.g. for g = 0.2, due to the repulsion of 
the bosons is observed. In this region the system is de¬ 
scribed by a superposition of lowest-band states which 
are predominantly of single-pair occupancy, e.g. |1, 2, l)j, 
|2, l,l)j, and double-pair occupancy, e.g. |2, 2,0) r For 
further increasing repulsion, e.g. g = 0.4, a trend towards 
the repopulation of the central well is noted. As we en¬ 
ter the strong interaction regime, e.g. g = 1.5, the state 
consists of a particle in the first excited-band being on 
a commensurate background of localized particles which 
lie in the zeroth band and the dominant ground state 
configuration is 11, 2, l)j, with I L = I R = (1,0,0) and 
Im = (1,1, 0). Finally, for strong interparticle repulsion, 
e.g. g = 3, the contribution from the higher-band states 
becomes more prominent and the corresponding ground 
state configuration is characterized by an admixture of 
zeroth- and excited-band states. 


III. Driven quantum dynamics 

This section is devoted to a detailed analysis of the 
bosonic dynamics in a driven optical lattice. At the be¬ 
ginning, a general overview of the effect of the driving on 
the finite bosonic ensemble with respect to the driving 
frequency is given. Subsequently, a Floquet analysis is 
employed in order to investigate the underlying single¬ 
particle physics. Finally, we focus on specific interaction 
effects. 


A. Dynamical response 

Let us explore the dynamical response or sensitivity 
of the system with respect to the driving frequency u>d- 
In order to investigate the stability of the system against 
the perturbations induced by the shaking [see Eq.(l)], we 
first analyse the fidelity [38j between the initial state and 
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FIG. 1: (a) Time evolution of the fidelity F U , D ( t ) as a function 
of the driving frequency ujd (measured in units of uir). (b) 
Time evolution of the expectation value of the Hamiltonian 
(2) (measured in units of the recoil energy Er ) for various 
driving frequencies uid = 0.4 (green thin dashed line), ujd = 
4.0 (black thick solid line), ujd = 4.5 (red thick dashed line), 
ujd = 5.25 (magenta thick dashed-dotted line), ujd = 11.0 
(blue thin dashed-dotted line), and u>d = 13.375 (light-blue 
thin solid line). The driving amplitude is A = 0.05, while the 
initial state corresponds to the ground state of four weakly 
interacting bosons with g = 0.1 confined in a triple-well. Time 
unit is UJft 1 . 


the state evolved at time t: F ulD {t) = | (41(0)1 4 '(t)) | 2 , 
where the dependence on ujd is implicit in the time 
evolved state 4>(t). Here we will consider a system of 
four bosons in a triple-well with g = 0.1, whose ground 
state (i.e., the initial state 41(0)) corresponds to a super- 
fluid state, as the filling factor is not commensurable and 
we do not encounter the formation of a Mott insulating 
state. In terms of its dominant spatial configuration our 
system initially consists (see also Sec.II.C) of two bosons 
in the middle well and two others each of them local¬ 
ized in one of the outer wells, i.e. the state 11, 2, l) l7 
with I L = Ir = (1, 0, 0) and Im = (2, 0, 0) has the most 
prominent contribution. Figure 1(a) shows F UD (t) as a 
function of the driving frequency ujd- The dynamics is 
characterized by three main regions with respect to ujd, 
where the system is driven far from the initial state, while 
for the remaining frequency regions (red sections in Fig. 
1(a)) the evolved state is essentially unperturbed by the 
driving. In the first region, between 4.0 < ujd < 5.5, 
the minimal overlap in the course of the dynamics drops 
down to 0.1, whereas in the second (7.0 < ud < 8.0) 
and third (10.0 < u>d < 15.0) regions the system maxi¬ 


mally departs from the initial state with a percentage on 
the order of 50% and 65%, respectively. The emergence 
of these dynamical regions strongly depends on the pa¬ 
rameters of the optical lattice. For instance, for smaller 
lattice depths the aforementioned regions will be wider, 
because of the smaller potential energy, which favors a 
possible deviation of the system from the initial state. 

Let us inspect the time evolution of the total energy 
E(t) = {${t)\H(t) |T(t)}. Figure 1(b) shows E{t) for 
various driving frequencies u>d ■ For driving frequencies 
where F UJD ~ 1 [e.g, ujd € {1,3}, see also Fig. 1(a)] 
the dependence of the energy on the driving frequency is 
weak and it is essentially constant during the time evolu¬ 
tion. On the other hand, for the regions where F U1D <C 1, 
E(t) increases initially and it shows an oscillatory be¬ 
haviour. In particular, for ujd = 4.5 the total energy 
exhibits an oscillatory (almost periodic) pattern which 
can also be observed in the corresponding fidelity evolu¬ 
tion. This driving frequency will be referred to in the fol¬ 
lowing as critical and denoted by uj^, that is, the driving 
frequency for which min tg roj’] F UD (t) is minimal. Indeed, 
as we shall see below, the most interesting dynamics of 
the system takes place close to this frequency. 

Finally, let us inspect the response of the sys¬ 
tem to the driving from a one-body perspec¬ 
tive via the single-particle density pi(x,t) = 
/ dx2...dxjv|4 , (x, X 2 ,..., xjv; t)| . Figure 2 illustrates 
the evolution of the one-body density for different 
driving frequencies wd, but with the same amplitude 
A. The driving leads to oscillations of the particles 
densities in every site. As it can be observed by having a 
glance at Fig. 2(a), the one-body density shows a weak 
response for driving frequencies away from the critical 
region ud £ [4,5.5], while for u>d = w C D [see Fig. 2(b)] 
we observe the periodic formation of enhanced density 
oscillations being accompanied by a broadening of each 
intra-well ensemble. The peculiar behaviour of the 
bosonic ensemble observed for ujd = is characterized 
by three processes and time scales: i) the internal 
fast oscillations of the density; ii) the large amplitude 
oscillations of the density in each well of period ~ 14; 
in) the tunneling between the wells with a period of 
about 200. All these features will be analyzed in detail 
in the following subsections both at the single particle 
and many-body level. 

B. Single particle dynamics 

Here we investigate to which extend the previously pre¬ 
sented results can be understood in the limit of zero inter¬ 
action among the particles by means of Floquet theory. 
Specifically, we are interested in two distinct features of 
the dynamics observed in Fig. 2(b): First, the on-site 
dynamics and, especially, its resonance-like dependence 
on the driving frequency wd, and second, the inter-well 
tunneling dynamics which is enhanced at certain values 
of u>d- 
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FIG. 2: Time evolution of the one-body density pi(x,t ) in a triple-well potential for different driving frequencies: (a) a in = 2.0 
(top panel) and (b) ojd = 4.5 (lower panel). The driving amplitude is fixed to the value A = 0.05, while the initial state 
corresponds to the ground state of four weakly interacting bosons with g = 0.1. The spatial extent of the lattice is expressed 
in units of fc^ 1 , while the time units are rescaled in terms of the driving period Tn- 


1. Floquet theory 

To be self-contained, we start by summing up the 
main notions of Floquet theory. Because of the temporal 
periodicity of the single particle Hamiltonian employed 
throughout this work [Eq. (2) with g = 0 and N = 1], 
every solution of the time-dependent Schrodinger equa¬ 
tion (TDSE) takes the form of a Floquet mode (FM) 
'F a (x,f) which in turn can be written as: ’F Q (a;,f) = 
e ~ie a t/hq> a (x,t) with the real quasi energy (QE) e a € 
[— hujp/2, +Tiuj d /2} and with <F a (x,f) = <F a (x,f + T D ) 
resp ecting the temporal periodicity of the Hamiltonian 
[39j . The FMs are eigenvectors of the time evolution op¬ 
erator over one driving period 

U{T d + t 0 ,t 0 )^ a {x,t 0 ) = (4) 

This property is of particular interest as it allows for a 
stroboscopic time evolution of an arbitrary initial state 
^(x,^) once the FMs of a system are know. To show 
this, we exploit the fact that the FMs constitute an or¬ 
thonormal basis for the solution space of the TDSE [40| 
and expand ^(x, to) at the initial time t = to as 

^(x,t 0 ) = ^2c a (to)^f a (x,t 0 ) (5) 

a 

with the corresponding coefficients C a (to). By applying 
the one period evolution operator U(Tp + t Ql t 0 ) on both 
sides of Eq. (5) for m times and by virtue of Eq. (4), we 
readily obtain the stroboscopic time evolution of T(x, to) 


as 

T(x, t 0 + mT) = J2 C a (t 0 )e- ie “ mTD / fi 4/ a (x, t 0 ). (6) 

Ot 

Numerically, we obtain the FMs for a given initial time 
to by calculating the eigenvectors of the one period evo¬ 
lution operator U(Tp+tn, to) [see Eq. (4)]. We refer the 
interested reader to ref. [4J| for a detailed description of 
the employed computational scheme. 

Finally, let us note that Eq. (6) already reveals some 
interesting features of the time evolution in periodically 
driven systems as we shall see in the following. Imag¬ 
ing that only a single FM, say 'F 0 (x, t), is populated. 
The stroboscopic evolution of the probability density is 
thus given as |\H(x, toTd)| 2 = |Co| 2 |'F Q ,(x,0)| 2 . Hence, 
|’F(x,f)| 2 is again periodic with period Tp and the only 
time dependence arises from the explicit time dependence 
of the FM, which is commonly referred to as ’’micro¬ 
motion” . This situation changes if the initial state popu¬ 
lates multiple FMs. In this case one encounters interfer¬ 
ence terms in |’F(x, mT£>)| 2 between the different FMs in 
the form of ~ e ™(e«-£/))TD/fi i Thus, the quasi energies 
e a in periodically driven systems play a comparable role 
in the time evolution as the energy eigenvalues do in time 
independent setups. 

2. On-site dynamics in the single well 

To begin with, we shall investigate the observed on¬ 
site dynamics [see Fig. 2(b)], and in particular their de¬ 
pendence on the driving frequency cop. To this end, we 
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FIG. 3: On-site dynamics for a single particle, (a) Popula¬ 
tions || 2 and |Ci| 2 of the two most populated FMs. (b) QE 
spectrum as a function of the driving frequency lod (measured 
in units of lor). Highlighted are the most (red) and second 
most (blue) populated FMs. The rectangular area indicates 
the narrow avoided crossings, while the circle highlights the 
area where a broad avoided crossing among the FMs appears 
with respect to the driving frequency, (c) In black is again 
the QE spectrum (same as in (b)). Additionally, we show 
the most (red) and second most (blue) populated states of 
the static, i.e. undriven, lattice. For comparison, we depict 
again the black rectangle at the same position as in (b). (d) 
Frequencies io a p of the on-site dynamics as a function of the 
driving frequency (see main text), (e) Same as (d), but in the 
triple well setup (shown is only the extract of small frequen¬ 
cies uJa/3 -C lod which corresponds to tunneling dynamics). In 
all panels A = 0.05. 


simplify the setup studied in Sec. III.A to just a single 
well of the lattice potential. Hence, the potential is given 
by V sh (x,t) = Vo sin 2 (x — Asin(wut)) for x £ (—7r,+7r] 
and we impose periodic boundary conditions at x = ±7r 
in order to mimic the situation in an extended lattice. 
We choose as initial state T(:r, 0) the single particle den¬ 
sity as shown in Fig. 2(b) at t = 0 within the central 
potential well. The time evolution is then obtained by 
expanding T(a:, 0) in terms of the FMs of the system and 
by making use of Eq. (6). As a result, we find that we 
can reproduce some of the main features of the on-site dy¬ 
namics shown in Figs. 2(a) and (b), namely, we observe 
resonantly enhanced on-site oscillations in an interval of 



FIG. 4: (a) Tunneling probability (see main text) I>;i I“ = 
|i (2,1, l|T(t)) | 2 with I L = (2, 0, 0) and I R = I M = (1, 0, 0) as 
a function of time for different driving frequencies lod = 0.5 
(blue solid line), lod = 4.5 (red dashed line) and lod = 11.0 
(black dashed-dotted line). The most significant contribu¬ 
tion of the interband tunneling mode is between the state 
|2,1, l)j [with II = (2, 0, 0) and Ir — Im = (1,0,0)] and the 
initial 11, 2, l)j [with Im = (2,0,0) and Ir — II = (1,0,0)]. 
(b) Inter-well tunneling probability I jv^ > ; i | 2 at resonance for 
different values of the interatomic interaction g = 0.1 (black 
dashed-dotted line), g = 0.5 (blue solid line) and g = 2.0 (red 
dashed line). In all panels A = 0.05. The time evolution is 
expressed in units of w^ 1 . 


the driving frequencies around ojd ~ 4.5. Following the 
discussion in ref. [42|] , further insight into this effect can 
be obtained by studying the population of the FMs by 
the initial state as a function of ujd- We therefore sort 
the FMs according to their overlap with the initial 
state and label the mode with the largest overlap as T 0 , 
the mode with the second largest overlap as Ti etc. In 
Fig. 3(a) the coefficients of the two most populated FMs, 

|Co | 2 and |Ci | 2 , are shown as a function of the driving fre¬ 
quency. Apparently, both at small frequencies (lod 4) 
and at large ones (ujd > 5.5) only a single FM is notably 
populated, while |Co | 2 and |Ci | 2 become comparable at 
distinct driving frequencies (e.g., at ujd ~ 5). According 
to our discussion above, in cases when |Co | 2 is close to 
one, and thus only a single FM is populated, the strobo¬ 
scopic time evolution, as given by Eq. (6), becomes, to a 
good approximation, time periodic with the period of the 
driving To- Note that this agrees with the observation 
of Fig. 2(a), that away from the resonance frequencies, 
the single particle density merely performs oscillations 
whose period matches Tp. This corresponds precisely to 
the previously described micro-motion arising from the 
explicit time dependence of the FM To (x,t). 

On resonance, when |Co| 2 ~ |Ci| 2 , the evolution of 
|T(x, 0)| 2 includes, besides the micro-motion, an inter¬ 
ference term between To and Ti, whose period is dic¬ 
tated by the corresponding quasi energies and is given 
by: T osc /To = htoo/i^i — eo)- Indeed, we find that this 
term is responsible for the observed on-site mode with a 
period of ~ 14 lattice oscillations [compare Fig. 2(b)]. 
Up to now, however, it is not yet clear why Ti is res- 
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FIG. 5: (a) Probability 173{> ; i 1 2 (see main text) for all 
bosons to be in the zeroth-band during the evolution for dif¬ 
ferent driving frequencies ton = 2.0 (blue solid line), ton = 4.5 
(black dashed line), ton = 10.25 (red dotted line) and cod = 
11.75 (green dashed-dotted line), (b) Comparison of differ¬ 
ent excitation scenarios at ujd = 4.375. The black dashed- 
dotted line refers to the probability for all the bosons to be in 
the zeroth-band while the blue dotted, red dashed, light-blue 
thick solid and magenta thin solid line refer to the probabil¬ 
ity to have one, two, three or four bosons, respectively, in 
the first-excited band, (c) Probability Ijv^ };i I^ f° r ah the 
bosons to be in the zeroth-band for ton = ton, but with dif¬ 
ferent interparticle repulsion g — 0 (black dashed-dotted), 
g = 0.1 (blue dotted line), g = 0.5 (red dashed-dotted line) 
and g = 2.0 (light-blue solid line). In all panels A = 0.05. 
The time evolution is expressed in units of to^ 1 . 


onantly populated at certain frequencies. In order to 
provide an answer to this question we follow the argu¬ 
mentation in ref. [42j and consider the dependence of 
the quasi energy spectrum on the driving frequency ton 
as shown in Fig. 3(b). Highlighted are the two most pop¬ 
ulated modes at each wjj (blue and red dots) revealing 
avoided crossings of these two modes at the frequencies 
where a resonant enhancement of \C \\ 2 was observed in 
Fig. 2(b). Hence, at these values of wd the FMs d>o and 
di i are resonantly coupled by the driving which results in 
an increase of |Ci| 2 and ultimately to the in Sec. III.A 
described on-site dynamics. 

In the following we provide insight into the question 
why we observe Floquet resonances at driving frequen¬ 
cies around lod ~ 4.5. Let us start by noting that, by 
means of appropriate unitary transformations, the single 
particle Hamiltonian with a potential as given in Eq. (1) 
can be recasted into the form: 

h 2 d 2 

H = + Lo sin 2 (fco®) + V D sin{co D t)x , (7) 

where the amplitude of the oscillating term is given by 
Vjj = mAto 2 D . That is the transformed Hamiltonian takes 
the form of a static lattice plus a time-dependent pertur¬ 
bation whose strength is determinded by Vd- For the 


used parameters of m = 1 and A = 0.05 and for the 
range of considered frequencies of 3 < wp < 6 we get 
that the amplitude Vd of the time-dependent term is of 
order one. Hence, it can be seen as a small perturbation 
compared to the static term of strength Vo = 15 and we 
can expect that the QEs of the driven lattice setup can 
be estimated by the actual energies of the undriven lat¬ 
tice. Resonances would than be expected whenever the 
energy difference between two notably populated eigen¬ 
states of the static system matches an integer multiple 
of Tilud- In fact we find that the energies of the three 
energetically lowest states of even parity are given by 
E 0 « 2.6, Ei ss 11.6 and E 2 = 16.9. Naivly, we expect 
driving induced resonances whenever the ground state is 
resonantly coupled to one of the excited states. Indeed 
we find E\ — Eq « 2 x 4.5 and E 2 — Eg ~ 3 x 4.8. Thus, 
following this line of arguments, at driving frequencies 
of approximately lod = 4.5 the ground state of the un¬ 
perturbed lattice is coupled via a 2 (3) photon process 
to the first (second) excited states. In order to justify 
this simplyfied picture we show the energies Eq and E\ 
on top of the QE spectrum of the driven lattice (see Fig. 
3(c)). Away from any resonances, both energies are al¬ 
most identical to the QEs of the corresponding Floquet 
states, so e.g. the red line is practically on top of an un¬ 
derlying black line. Closer to the resonance region, we 
see of course deviations of the QEs from the mere en¬ 
ergies of the undriven lattice as the different states are 
coupled by the driving. 

Finally, Fig. 3(c) provides an overview over the possi¬ 
bly observed frequencies in the on-site dynamics at vari¬ 
ous driving frequencies. Shown are the frequencies asso¬ 
ciated to all possible interference terms between the FMs 
weighted by their overlap with the initial state. More 
precisely, we calculate tu a 0 = ( e a — ep)/Ti for all pairs 
of FMs at a given driving frequency and determine the 
colour coding by computing the product \C*Cg\. Hence, 
the frequency tu a p appears in Fig. 3(c) only when both of 
the corresponding FMs ’F Q and 4^ have appreciable over¬ 
lap with the initial state. In agreement with the discus¬ 
sion concerning Fig. 3(b) we observe pronounced on-site 
oscillations only within an interval of driving frequencies 
4.0 u>d 5.5. In particular, the two narrow avoided 
crossings around to c D « 4.5 [see the rectangular in Fig. 
3(b)] yield a low frequency on-site dynamics, whereas 
the comparably broad avoided crossing at lod ~ 5 [see 
the circle in Fig. 3(b)] results in a much faster on-site 
oscillation. 


3. Tunneling dynamics in the triple well 

Besides the on-site dynamics, Fig. 2(b) revealed a 
pronounced tunneling between the lattice sites at cer¬ 
tain driving frequencies. Similar to the previous sec¬ 
tion, we analyze this effect in the following by apply¬ 
ing Floquet theory for the single particle dynamics. We 
choose the same setup as before, that is, V s h{x,t ) = 















Vosin 2 (x — Asiniujot)), with the same initial state, i.e. 
essentially a Gaussian centered around the potential well 
at x = 0, but with the difference that the periodic bound¬ 
ary conditions are imposed at x = ±27r (instead of at 
x = ±7r as we did before). In this way we allow for 
tunneling of the wave packet into the two neighbouring 
lattice sites. As for the on-site dynamics, we provide an 
overview over the observable frequencies in the temporal 
evolution in Fig. 3(d) [in close analogy to Fig. 3(c)]. 
Note that, since the tunneling dynamics observed in Sec. 
III.A occurs on much longer timescales as compared to 
the on-site dynamics, we only show the extract of the 
regime of small frequencies, i.e. -C uid- Further¬ 
more, because no on-site dynamics occurs with timescales 
matching the extremely small frequencies of < 0.02, all 
the frequencies depicted in Fig. 3(d) are indeed associ¬ 
ated with an inter-wcll tunneling mode. In accordance 
with the observation made in the many-particle simu¬ 
lations [cf. Fig. 2(b)] we observe a strong increase of 
the frequencies associated with the tunneling dynamics 
in the range of driving frequencies of 4 < lod 1$ 5.5. 
Away from this resonance, for example at ojd = 2.5, 
the only notable tunneling mode corresponds to an inter¬ 
ference term of two FMs which oscillates with a period 
of T osc /T d ~ 3300 and could therefore not be observed 
in the simulations performed in Sec. III.A. Within the 
regime of resonant driving, e.g. at u>d = 4.5, the fre¬ 
quency of the tunneling mode is increased strongly and 
the associated oscillation period becomes T osc /Td « 200 
matching the observed tunneling mode in the weakly in¬ 
teracting regime [cf. Fig. 2(b)]. 

C. Interband tunneling and excitation processes 

In the previous section we have shown that most of the 
features of the (effective) single-particle dynamics of Fig. 
2 can be explained via a non-interacting Floquet theory. 
As we shall see now, however, the full dynamics presents a 
rich excitation spectrum ascribable to the particles inter¬ 
action, especially in the strong interaction regime. Thus, 
we investigate the tunneling and excitation probabilities 
of the dominant particle configurations, for different driv¬ 
ing frequencies wp, by means of the multiband expansion 
introduced in Sec. II.C. More precisely, we compute and 
analyze the probabilities, during the dynamics, defined 
as 

iq7v i}: i| 2 = |i(^i,^ 2 ,Ar 3 |vl/(t))| 2 . (8) 

(k) 

The case Im = 0 V k > 1 refers to the lowest-band 
inter-well tunneling dynamics. The initial state of the 
system corresponds to the ground state of four weakly 
interacting bosons with g = 0.1 in a triple well, while the 
dominant number state configuration (see also Sec.II.C) 
is 11,2, l)j with I L = Ir = (1,0,0) and I M = (2,0,0). 
In this way, a lowest-band tunneling process can take 
place among the initial state and: a) another state of 
single-pair occupancy, e.g. |2,1, l) z (II = (2,0,0) and 


Im = Ir = (1,0,0)), b) a state with double-pair oc¬ 
cupancy, e.g. 12,2,0^ ( I L = I M = (2,0,0) and I R = 
(0,0,0)), c) a state with triple occupancy, e.g. |3,1, 0) z 
(II = (3,0,0), I M = (1,0,0), I R = (0,0,0)) or d) a 
state with quartic occupancy, e.g. |4, 0,0)j (II = (4,0, 0) 
and Im = Ir = (0,0,0)). However, the from the system 
prefered tunneling processes form a hierarchy according 
to the energetical difference between the initial and final 
state. For instance, a tunneling process to another state 
of single-pair occupancy will be more preferable than to a 
state of double-pair occupancy etc. Figure 4(a) shows the 
tunneling probability to the energetically closest number 
state, which is |2,l,l)j (or 11,1,2) x ) with II = (2,0,0) 
and I M = I R = (1,0,0) (or I L = I M = (1,0,0) and 
Ir = (2,0,0)), i.e. |T»{iv i };i| 2 = |i < 2 , i, l| 1 T'(it)> | 2 with 
II = (2,0, 0) and Im = Ir = (1,0, 0), for various driving 
frequencies. As it is shown, for lod < w c D this tunnel¬ 
ing mode has a small amplitude and it is quite insensi¬ 
tive to wn as intuitively expected from the fact that the 
evolved-state is essentially unperturbed by the driving 
(see also Fig. 1(a)). For wd ~ ocfj, however, the ampli¬ 
tude of the oscillations is significantly larger indicating 
an enhancement of the tunneling (see also Figs. 2(b) and 
3(d)), whereas when F UD <C 1 and u>d > (wo = 11.0 
curve in Fig. 4(a)) the oscillations occur with a larger 
period. The fact that it oscillates with a larger time pe¬ 
riod can be traced back to the behaviour of the fidelity 
at wn = 11. Indeed, for short times (ojd = 11.0) the 
system stays in the initial ground state and after some 
time the fidelity starts to decrease, differently from the 
situation at ujd = 4.5, where the system deviates from 
the initial state on much shorter time scales. Concerning 
the remaining tunneling modes, i.e. tunneling to higher 
energetical states that belong to the lowest-band (see dis¬ 
cussion above), they are negligible as they provide a very 
small contribution even for u>d = The latter has al¬ 
ready been seen in the last subsection, but it can also be 
shown with the use of the multi-band analysis. On the 
other hand, Fig. 4(b) presents again the tunneling prob¬ 
ability |H{iv i } : i| 2 for the energetically closest lowest-band 
states (i.e. same as Fig. 4(a)) when the driving frequency 
is at resonance for different interaction strengths g. For 
weak to intermediate interactions the tunneling ampli¬ 
tude decreases and for strong interactions, e.g. g = 2.0, 
a destruction of the tunneling is observed for long time 
scales. 

Now, let us consider the excitation dynamics. In this 
(k) 

case it holds Im 0 for k > 1. To this aim, we have an¬ 
alyzed the probability of finding all the four bosons in the 
zeroth-band. The latter, can be expressed via eq.(7) as 
\B{n & ,i \ 2 = Ei\i(Ni,N 2 ,N 3 \*(t))f = EilCwbil 2 , 
where the summation is performed over the excitation 
indices I = (/l, Im,Ir) which, in terms of the multiband 
expansion, obey the constraints 1 ^ +1$ +I R ' > = N and 
1^ = I$ = I= 0 for all j > 1. In particular, Fig. 
5(a) shows the probability \B{ N .yj \ 2 for all the bosons 
to reside in the zeroth-band for various driving frequen- 


9 


cies u>d and a fixed amplitude A = 0.05 during the time 
evolution. At the critical driving frequency a complete 
depopulation of the zeroth-band at some specific time 
intervals is observed. In particular, this probability ex¬ 
hibits revivals, which are connected with the enhance¬ 
ment of the (amplitude) oscillations of the single-particle 
density [see also Fig. 2(b)]. On the other hand, for driv¬ 
ing frequencies different from the critical frequency the 
respective probability for all the bosons to occupy the 
zeroth-band is rather large and is indeed dominant. How¬ 
ever contributions from excited configurations cannot be 
neglected, especially in the regions 7.0 < ojd < 8.0 and 
10.0 < ojd < 15.0, where the system significantly de¬ 
parts from the initial state [see also Fig. 1(a) and Fig. 
5(a) red dashed line]. Furthermore, Fig. 5(b) presents 
the probability, at the critical driving frequency, to ob¬ 
tain a state of No < 4 particles in the first-excited band 
and the remaining to be in the zeroth-band. The latter 
can be expressed as |Q{jy 4 }:il 2 = Si |C{A5} ; i| 2 where the 
summation index I = (II, ImJr) obeys the constraints 


r(i) 


+ Im + Ir =N-N 0 „ 


r(b 


r(2) t ( 2) 

1 L ~~ 1 M 


= 4 2) = 


and l]p + 1$ + I = 0 for all j > 2. Indeed, the 
interplay between the four possible excitation scenarios 
from the zeroth to the first excited-band (i.e. one-particle 
excitation, two-particle excitation etc) in the course of 
the dynamics is illustrated in a transparent way. It is 
observed that the complete depopulation of the zeroth- 
band is mainly accompanied by the excitation of three 
or all the four bosons in the first-excited band. For long 
evolution times the zeroth-band possesses a low popu¬ 
lation and states with one or two bosons in the first 
excited-band are mainly populated. The states with the 
most significant contribution are of the type 11,2, l)j with 
II = Ir = (0,1,0) and I M = (0, 2, 0) or I M = (1,1,0). 
We note that a small contribution comes from the state 
11,2, l)j with I L = I R = (0,1,0) and I M = (0,1,1). 
This clearly shows that the most prominent excitation 
process in our system originates from the energy differ¬ 
ence between each of the above states and |1, 2,1) 7 with 
II = Ir = (1,0,0) and Im = (2,0,0), namely the (ini¬ 
tial) ground state configuration. 

Finally, in order to explore the impact of the interac¬ 
tions on the dynamics, Fig. 5(c) shows the probability 
\ B {Ni}:l\ 2 for long evolution times for all the bosons to 
be in the zeroth-band for different interparticle repul¬ 
sion at the driving frequency lod = For the non¬ 
interacting case the population of the zeroth-band shows 
revivals even for long time scales, while, as the interac¬ 
tion strength is turned on, the corresponding probability 
presents a decaying envelope. This envelope behaviour 
is a pure effect of the interactions and reflects also the 
initial ground state configuration (see the discussion in 
Sec. II.C) which strongly depends on the interparticle 
interactions. As it can be seen for increasing repulsion 
between the particles the probability for the system to 
remain in the zeroth-band, in the course of the dynam¬ 
ics, decays on increasingly shorter time scales and the 
system is dominated by different types of excitations, as 


expected intuitively. 
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FIG. 6: Local one-body density spectrum pl{ to) (for the left 
well) as a function of the driving frequency ujd (measured 
in units of lor). The driving amplitude has been chosen 
A = 0.05. Inset: the spectrum of the intra-well oscillations 
calculated via Apn(t) (see also main text). 



FIG. 7: (a) Profile of the resonance for various driving am¬ 
plitudes A= 0.01 (blue solid line), A= 0.05 (black dashed- 
dotted line) and A= 0.1 (red dashed line) obtained from the 
m i n te[o,T] Po(t) and T being some fixed long evolution time, 
as a function of the driving frequency, (b) Same as (a) with 
A = 0.05, but for different barrier heights Vo= 9.0 (red solid 
line) and Vo= 12.0 (blue dashed line). The system consists of 
four bosons confined in a triple-well with interparticle inter¬ 
action g = 0.1. 


D. Characteristics of the resonant behaviour 

To characterize the overall process with respect to the 
driving frequency, we compute the spectrum of the local 
one-body density 


Pa(w) = - f dtp a (t)e zut , (9) 

* J o 

where p a (t) = J^ a dxp\{x,t ) denotes the spatially over a 
single well integrated single-particle density at every time 
instant t. The index a = L, M, R corresponds to the left, 
middle or right well respectively, whereas the limits of 
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FIG. 8: Momentum distribution of the one-body density as a 
function of time (measured in units of uq) 1 ) for g = 0.1 and 
different driving frequencies (a) before the critical frequency 
lod = 2.0, (b) at the critical frequency ujd = uib = 4.5 an< l 
(c) at ujd = 8.0. (d) The case of strong interparticle repulsion 
for g — 2.0 and u>d = ub- The horizontal axis represents the 
lattice momenta in units of the inverse lattice vector ko = n/l. 
In all panels A = 0.05. 


the wells are denoted by d a , d' a . Note that in the present 
case all the components of p a {uj), i.e. pl{ui), Pm{w) and 
Pr(uj), are equivalent due to the considered large lattice 
depths and the employed driving scheme which enforces 
the bosons among different wells to oscillate in-phase. 
Figure 6 shows the above spectrum, where five dominant 
branches (denoted as (1) to (5) in the Figure) can be ob¬ 
served. The lowest branch denoted as (1) in Fig.6 (in the 
range w € [0,0.02]) refers to the intraband tunneling be¬ 
ing restricted to the energetically closest number states 
e.g. from 11,2,1)! (J M = (2,0,0), I L = Ir = (1,0,0)) 
to 12,1,1)! (I L = (2,0,0), J M = Ir = (1,0,0)). This 
branch is hardly visible in Fig.6 due to the presented wide 
range of frequencies that have been taken into account in 
order to visualize all the dynamical frequencies of the 
system. In addition, the next lowest branch (denoted as 
(2)) at ujd € [4,5] and w € [0.05,1] corresponds to the 
large amplitude density oscillations [see also Fig. 2(b)]. 
These mode frequencies have been already predicted via 
the Floquet analysis in Sec. III.B [see Fig. 3(c) and (d)]. 
To investigate in some detail the intra-well wavepacket 
dynamics the quantity A p a (t) = p a p(t) — p a ,i{t) is 
employed. Here, each well is divided from the cen¬ 
ter into two equal parts, namely left and right, with 
Pa,i(t), Pa, 2 (t) being the corresponding integrated den¬ 
sities at time t. The index a = L,M,R stands for the 
left, middle and right well, respectively. To determine 
the frequencies of this mode we calculate the spectrum 
A pl(w) = i / dt/S.pL(t)e lult . The inset of Fig. 6 presents 
the corresponding spectrum, thus showing the emergent 
frequencies of the intra-well oscillations as a function of 
the driving frequency ujd • We observe that the spectrum 
A pl(w) follows the evolution of the upper three branches 


(denoted by (3), (4) and (5)) of the spectrum of pl{w), 
whereas in the region of the resonance the intra-well os¬ 
cillation measured via A pl(Jj) features a beating dynam¬ 
ics, as expected. Hence, away from the region around 
the critical driving frequency the generated dipole mode 
possesses three different frequencies, while close to uj c d 
the intra-well dynamics come into a resonance. There¬ 
fore, one can induce this resonant intra-well dynamics by 
adjusting the driving frequency. Finally, let us comment 
on the existence of some higher frequency components, 
e.g. branch (6) in Figure 6, which correspond to very 
fast intrawell oscillations (i.e. uj <C wo) and possess a 
low amplitude (in comparison to the previous branches 

(l)-(5))- 

In turn, we shall visualize the above mentioned reso¬ 
nance and inspect how it depends on the lattice param¬ 
eters. To this aim, the minimal occupancy, during the 
evolution time T, of the zeroth-band min tg [ 0 ,T] Po(t) = 
min te[o,T] Ei |l (Ni,N 2 ,N 3 \^!(t)) | 2 , with the energetical 
indices 4 ° + 7 $ + 7 £> = N and I ( L j) = I $ = 1$ = 0 
for every j > 1 is used. Employing the above quantity 
one can show that far from resonance there are regions 
with non-negligible excitations i.e. min tg [ 0i T] Po < 1 
[e.g. at u jd = 11-0 see also Fig. 1(a)] as well as re¬ 
gions where min tg r 0i T] Po ~ 1 [e.g. <^d = 2.0 in Fig. 
1(a)]. Now let us analyze the dependence of min te ro ;T i P 0 
on the driving frequency around Firstly we study 
the dependence of the resonance on the driving ampli¬ 
tude. In Figure 7(a) we show for an increasing driv¬ 
ing amplitude the minimum of min te [o ;T i P 0 as a func¬ 
tion of the frequency u>d which broadens and eventu¬ 
ally reaches zero, meaning that the zeroth-band has been 
complete depopulated [see also Fig. 5(a)]. On the other 
hand, for small amplitudes the value of the minimum of 
min t6 r 0jT ] Po is non-zero and in the limit A —> 0 its de¬ 
pendence on the driving frequency disappears. Instead, 
in Fig. 7(b) we show how the minimal population of the 
zeroth-band (rriin tg [o ; T] Po) varies as a function of the lat¬ 
tice depth. For an increasing lattice depth it is known 
that the energy gaps among the different energy levels 
become larger. This phenomenon can intuitively be un¬ 
derstood in terms of a tight-binding approximation. For 
simplicity let us assume only a nearest neighbour cou¬ 
pling J oc f dxW s {x)[ ^ + Vo sin 2 (a;)]IF s+ i(a;) between 
the sites s and s + 1, where W s (x) are the on-site lo¬ 
calized Wannier states. Then, within this approxima¬ 
tion, which is valid for a relatively deep potential, the 
resulting eigenvalues are E k -i = E° n ~ stte - 2 J cos(4fj) 
{k = 1, 2,..., N), where E° n ~ slte are the on-site energies. 
Thus, the resonance can be tuned at will, i.e. for a de¬ 
creasing lattice depth the is negatively shifted, as it 
is confirmed by the numerical results of Fig. 7(b). Fi¬ 
nally, let us comment on the dependence of the position 
of the resonance on the interparticle interaction strength 
g. Indeed, in order to investigate whether there is such 
a dependence various interaction strengths (for the same 
particle number N = 4), e.g. g = 0.1, g = 1.0 and 
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FIG. 9: Time evolution of the one-body density pi(x,t ) in a twelve-well potential for c od = 4.5. The driving amplitude is fixed 
to the value A = 0.05, while the initial state corresponds to the ground state of five weakly interacting bosons with g = 0.1. 
The spatial extent of the lattice is expressed in units of fcf, while the time units are rescaled in terms of the driving period 


T D . 


g = 3.0, have been considered (omitted here for brevity) 
and it was found that the position of the resonance is 
essentially unaffected. 

In the following, let us inspect the momentum distri¬ 
bution with varying driving frequency with the aim of 
understanding whether signatures of a parametric am¬ 
plification of matter-waves can be observed. The mo¬ 
mentum distribution is a routinely employed observable 
in atomic quantum gases experiments as it is accessible 
via time-of-flight measurements Q ■ This quantity can be 
calculated as the Fourier transformation of the one-body 
reduced density matrix as 

n{k, t) = j J dxdx'p\(x, x'\t)e~ lk ^ x ~ x \ (10) 

Here p\{x,x'\t) denotes the one-body reduced density 
matrix, being obtained by tracing out all the bosons but 
one in the density of the iV-body system. The panels 
8(a)-(c) of Fig. 8 present the time evolution of the mo¬ 
mentum distribution for different driving frequencies be¬ 
fore, on, and after the resonance. As it can be noted, 
exactly at the resonance the momentum distribution ex¬ 
hibits a special pattern, that is, some additional lattice 
momenta are periodically activated during the dynamics. 
In particular, it is observed that the modes ±4p ~ ±1.57, 
±k 0 ~ ±3.14, ±^f ~ ±4.713 are populated, whereas out 
of resonance only the ±fco modes are significantly pop¬ 
ulated. The population of the ±fco/2, ±3fco/2 modes at 
ujd = is reminiscent of the parametric amplification 
of matter-wave phenomenon, as observed experimentally 
in ref. pi} . However, an exact correspondence with ref. 
[2II ] cannot be made due to the very different setup of our 
system, i.e. its finite size and the hard wall boundaries. 
A detailed study of this process, also for higher parti¬ 
cle numbers and lattice potentials, would be desirable, 
but it is clearly beyond the scope of this work. Further¬ 
more, Fig. 8(d) shows the momentum distribution at res¬ 
onance, but for a strong interparticle repulsion g = 2.0. 
The expected periodic pattern for large evolution times 


is blurred as an effect of the strong interaction which 
decreases the degree of coherence. 

Finally, in order to demonstrate that our findings are 
of general character we investigate a larger lattice sys¬ 
tem with a filling factor smaller than unity. Specifically, 
the case of five bosons in a twelve-well finite lattice has 
been considered. Concerning the ground state with fill¬ 
ing factor v < 1, the most important aspect is the spatial 
redistribution of the atoms as the interaction strength in¬ 
creases. Indeed, as the repulsion increases from the non¬ 
interacting to the weak interaction regime the atoms are 
pushed from the central to the outer sites which gain and 
lose population in the course of increasing g. 

In the following, the shaking dynamics applied at 
t = 0 to the ground state of the five bosons which are 
trapped in the twelve-well potential in the weak interac¬ 
tion regime (g = 0.1) is explored. The emergent non¬ 
equilibrium behaviour shows similar characteristics as in 
the previous setup with filling v > 1, i.e. the occurrence 
of an intrawell dipole and an interwell tunneling mode. 
Interestingly, at the same frequency u>d = = 4.5 a 

resonance of the intra-well dynamics is observed. Fig¬ 
ure 9 presents the one-body density evolution exactly at 
the critical point ujd- As in the case for setups with 
filling v > 1, the formation of enhanced density oscilla¬ 
tions at each site is observed, which is in relation to the 
time periods where the zeroth-band is completely depop¬ 
ulated during the evolution. Employing a corresponding 
number state analysis the significant contribution of two 
kinds of number states has been confirmed: a) Either 
if ± ... ± if = N - 1, if = ... = if = 0 and one 
with if = 1 for k = 1,..., 12 or b) if ± ...±if = N— 1, 
if = ... = if = 0 and a certain if = 1 for 
k = 1,..., 12. Notice that the same kind of number states 
have been found to contribute significantly also in the 
dynamics of four bosons in the triple-well. The above 
mentioned observations suggest a generalization of the 
observed phenomena to larger systems as well. Indeed, 
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the same shaken scheme has been tested in different sys¬ 
tems (omitted here for brevity), e.g. 10 bosons in a triple¬ 
well, 6 bosons in five wells etc, confirming that the above 
observed resonant-like behaviour of the bosonic ensemble 
occurs in each setup. 

IV. Conclusions and Outlook 

The correlated non-equilibrium quantum dynamics of 
few-body bosonic ensembles induced by the driving of 
a finite-size optical lattice has been investigated. Our 
work focuses particularly on the regimes of large lattice 
depths and small driving amplitudes. This choice has 
been made in order to limit the degree of excitations that 
would otherwise lead to heating processes. Starting from 
the ground state of a weak or strongly interacting small 
ensemble, we have examined in detail the time evolution 
of the system induced by periodically driving the opti¬ 
cal lattice. We find that the dynamical evolution of the 
system is governed by two main modes: the inter-well 
tunneling and the intra-well dipole-like mode. The dy¬ 
namical behaviour of the system in the non-interacting 
regime has been firstly analyzed via Floquet theory, that 
is, at the single-particle level, providing an accurate in¬ 
terpretation of the observed processes. For large particle 
numbers and large interaction strengths, however, such 
a single-particle description was not anymore sufficient 
to provide an exhaustive explanation of the observed dy¬ 
namics, and a multi-band Wannier number state expan¬ 
sion has been employed. 

The inter-well tunneling mode is weak as a consequence 
of the deep optical lattice and the small driving ampli¬ 
tude. On the other hand, the local dipole mode has 
been identified from the intra-well oscillations of bosons 
in the individual wells. Remarkably enough, it has been 
found that by tuning the driving frequency the intra-well 
dynamics experiences a resonant-like behaviour. This 
is manifested e.g. by the enhanced oscillations in the 
one-body density evolution or from the periodic popu¬ 
lation of additional lattice momenta in the momentum 
distribution of the one-body density. Additionally, on 
a single-particle level in terms of Floquet theory, it has 
been shown that in the proximity of the resonance the 
first two FMs possess the main contribution, while away 
from resonance the dynamics can be described with the 
inclusion of the first FM. To explain the enhanced pop¬ 
ulation of the second FM at resonance the correspond¬ 
ing quasienergy spectrum has been employed, revealing 
avoided-crossings between the first two FMs at certain 
driving frequencies. To obtain the frequencies which refer 
to the on-site and tunneling dynamics, the correspond¬ 
ing frequencies associated with the interference terms be¬ 
tween the FMs have been employed showing pronounced 
on-site oscillations and an enhancement of the inter-well 
tunneling mode in the vicinity of the resonance. Consid¬ 
ering an ensemble of few-bosons we examined the influ¬ 
ence of the interatomic interactions both for the inter- 


and intra-well generated modes. Indeed, it has been 
found that the repulsion affects each of the aforemen¬ 
tioned modes, yielding a destruction of the inter-well tun¬ 
neling for strong interactions and an enhancement of the 
excitations (i.e. the contribution of higher-band states). 
Moreover, in the spectrum of the local one-body density 
with respect to the driving frequency all the relative dy¬ 
namical frequencies, e.g. on-site oscillations and tunnel¬ 
ing period have been identified. Finally, the occurrence 
of the above resonance seems to be universal in a peri¬ 
odically driven lattice as it is independent of the filling 
factor, the boundary conditions or the interparticle re¬ 
pulsion. 

We would like to underline that, contrarily to related 
studies based e.g. on effective model Hamiltonians or lat¬ 
tice calculations with tensor network methods, our many- 
body analysis based on the ab initio MCTDHB method 
has the advantage to provide the complete system wave- 
function in space and time. Thus, it enables us to ac¬ 
curately identify the involved intra- and inter-well band 
excitations. 

Let us comment on possible future investigations. Al¬ 
though in the present work we did not employ the multi¬ 
layer structure of the ML-MCTDHB method, our ab- 
initio approach is well suited to describe the dynamics of 
multi bosonic species. Given this, a first natural exten¬ 
sion would be to study the driven dynamics of mixtures 
consisting of different bosonic species in order to unravel 
the induced excitation modes or to device schemes for se¬ 
lective transport of an individual bosonic component. In 
relation to the present study, it would be interesting to 
simulate the parametrical amplification of matter-waves 
with interesting applications, like the generation of four- 
wave mixing, entanglement production, but also for fun¬ 
damental tests of quantum mechanics with massive par¬ 
ticles like the Hong-Ou-Mandel experiment, as recently 
performed with a Bose-Einstein condensate [23| . 

A. Appendix A: The Computational Method: 

ML-MCTDHB and MCTDHB 

Our computational approach to solve the many-body 
Schrodinger equation of the interacting bosons relies 
on the Multi-Layer MultiConfiguration Time-De pen dent 
Hartree method for Bosons (ML-MCTDHB) [43l . |44| 
which constitutes an ab-initio method for the calcula¬ 
tion of stationary properties and in particular the non¬ 
equilibrium quantum dynamics of bosonic systems of dif¬ 
ferent species. For a single species it reduces to MCT¬ 
DHB which has been established in refs. [32, HH, |45[ 
and applied extensively [45l-l48j. The wavefunction is 
represented by a set of variationally optimized time- 
dependent orbitals which implies an optimal truncation 
of the Hilbert space by employing a time-dependent mov¬ 
ing basis where the system can be instantaneously opti¬ 
mally represented by the corresponding time-dependent 
permanents. To be self contained let us briefly introduce 
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the basic concepts of the method and discuss the main 
ingredients of our implementation. 

Within the MCTDHB method the time-dependent 
Schrodinger equation {iTid t — H) ^(x,t) = 0 is solved 
as an initial value problem |H/(0)} = I'I'o)- The many- 
body wavefunction which is expanded in terms of the 
bosonic number states |ni,n 2 , based on time- 

dependent single-particle functions (SPFs) | (j>i(t)), i = 
1, 2, M, reads 

i*«)> = E Cn(t) |ni,n 2 , (Al) 

n 

Here M is the number of SPFs and the summation n is 
over all the possible particle combinations rii such that 
the total number of bosons is conserved and equal to N. 
To determine the time-dependent wave function l’P(t)) 
we need the equations of motion for the coefficients (t) 
and of the SPFs | <f>i(t)). Following the Dirac-Frenkel 
[itl [Ho] variational principle i.e. (d^\idt — H l^) = 0 
we end up with the well-known MCTDHB equations of 
motion [32, [3J, [45], [5l| consisting of a set of M non-linear 
integrodifferential equations of motion for the orbitals 
which are coupled to the li near equations of 

motion for the coefficients. 

For our numerical implementation a discrete variable 
representation (DVR) for the SPFs and a sin-DVR, 
which intrinsically introduces hard-wall boundaries at 
both edges of the potential, has been employed. The 


preparation of the initial state has been performed by us¬ 
ing the so-called relaxation method in terms of which one 
obtains the lowest eigenstates of the corresponding 771 - 
well setup. The key idea is to propagate some trial wave 
function by the non-unitary operator e~ Hr . This 

is equivalent to an imaginary time propagation and for 
r —>• 00 , the propagation converges to the ground state, 
as all other contributions (i.e., e ~ EnT ) are exponentially 
suppressed. In turn, we periodically drive the optical lat¬ 
tice and study the evolution of ’f , (xi, x 2 , Xn\ t) in the 
m-well potential within MCTDHB. To ensure the con¬ 
vergence of our simulations we have used up to 9 single 
particle functions thereby observing a systematic conver¬ 
gence of our results for sufficiently large spatial grids. 
An additional criterion that confirms the achieved con¬ 
vergence is the population of the lowest occupied natural 
orbital kept in each case below 0.1%. 

Acknowledgments 

S.M. thanks the Hamburgisches Gesetz zur Forderung 
des wissenschaftlichen und kiinstlerischen Nachwuchses 
(HmbNFG) for a PhD Scholarship. P.S gratefully ac¬ 
knowledges funding by the Deutsche Forschungsgemein- 
schaft (DFG) in the framework of the SFB 925 ’’Light 
induced dynamics and control of correlated quantum sys¬ 
tems”. A.N. gratefully acknowledges discussions with 
Klaus Mplmer related to four-wave mixing. 


[1] N. Goldman, and J. Dalibard, Phys. Rev. X 4, 031027 
(2014). 

[2] N. Goldman, J. Dalibard, M. Aidelsburger, and N. R. 
Cooper, Phys. Rev. A 91 , 033632 (2015). 

[3] O. Morsch, and M. Oberthaler, Rev. Mod. Phys. 78 , 179 
(2006). 

[41 I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 
80 , 885 (2008). 

[5] M. Olshanii, Phys. Rev. Lett. 81 , 938 (1998). 

[6] R. Grimm, M. Weidemuller, and Y. B. Ovchinnikov, Adv. 
At. Mol. Opt. Phys. 42 , 95-170 (2000). 

[7] L. Santos, M. A. Baranov, J. I. Cirac, H. U. Everts, 
H. Fehrmann, and M. Lewenstein, Phys. Rev. Lett. 93, 
030601 (2004). 

[8] S. Inouye, J. Goldwin, M. L. Olsen, C. Ticknor, J. L. 
Bohn, and D. S. Jin, Phys. Rev. Lett. 93, 183201 (2004). 

[9] T. Kohler, K. Goral, and P. S. Julienne, Rev. Mod. Phys. 
78, 1311 (2006). 

[10] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. 
Mod. Phys. 82 , 1225 (2010). 

[11] M. Lewenstein, A. Sanpera, and V. Ahufinger, Ultracold 
Atoms in Optical Lattices: Simulating quantum many- 
body systems, (Oxford University Press, 2012). 

[12] D. I. Choi, and Q. Niu, Phys. Rev. Lett. 82 , 2022 (1999). 

[13] M. B. Dalian, E. Peik, J. Reichel, Y. Castin, and C. 
Salomon, Phys. Rev. Lett. 76 , 4508 (1996). 

[14] O. Morsch, J. H. Miiller, M. Cristiani, D. Ciampini, and 
E. Arimondo, Phys. Rev. Lett. 87, 140402 (2001). 


[15] E. Peik, M. B. Dahan, I. Bouchoule, Y. Castin, and C. 
Salomon, Phys. Rev. A 55, 2989 (1997). 

[16] M. Cristiani, O. Morsch, J. H. Muller, D. Ciampini, and 
E. Arimondo, Phys. Rev. A 65, 063612 (2002). 

[17] S. R. Wilkinson, C. F. Bharucha, K. W. Madison, Q. Niu, 
and M. G. Raizen, Phys. Rev. Lett. 76 , 4512 (1996). 

[18] Q. Niu, X. G. Zhao, G. A. Georgakis, and M. G. Raizen, 
Phys. Rev. Lett. 76 , 4504 (1996). 

[19] C. Sias, H. Lignier, Y. P. Singh, A. Zenesini, D. Ciampini, 
O. Morsch, and E. Arimondo, Phys. Rev. Lett. 100 , 
040404 (2008). 

[20] A. Eckardt, C. Weiss, and M. Holthaus, Phys. Rev. Lett. 
95 , 260404 (2005). 

[21] N. Gemelke, E. Sarajlic, Y. Bidel, S. Hong, and S. Chu, 
Phys. Rev. Lett. 95, 170404 (2005). 

[22] K. M. Hilligspe, and K. Mplmer, Phys. Rev. A 71, 041602 
(2005). 

[23] R. Lopes, A. Imanaliev, A. Aspect, M. Cheneau, D. Bo- 
iron, and C. I. Westbrook, Nature 520 (2015). 

[24] W. Zheng, and H. Zhai, Phys. Rev. A 89, 061603 (2014). 

[25] H. Lignier, C. Sias, D. Ciampini, Y. Singh, A. Zen¬ 
esini, O. Morsch, and E. Arimondo, Phys. Rev. Lett. 99, 
220403 (2007). 

[26] J. Struck, C. Olschlager, M. Weinberg, P. Hauke, J. Si- 
monet, A. Eckardt, M. Lewenstein, K. Sengstock and P. 
Windpassinger, Phys. Rev. Lett. 108 , 225304 (2012). 

[27] C. V. Parker, L. C. Ha, and C. Chin, Nature Phys. 9, 




14 


769-774 (2013). 

[28] S. Choudhury, and E. J. Mueller, Phys. Rev. A 90 , 
013621 (2014). 

[29] C. Strater, and A. Eckardt, Phys. Rev. A 91 , 053602 
(2015). 

[30] A. Iucci, M. A. Cazalilla, A. F. Ho, and T. Giamarchi, 
Phys. Rev. A 73 , 041608 (2006). 

[31] P. I. Schneider, and A. Saenz, Phys. Rev. A 85 , 050304 

( 2012 ). 

[32] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, J. 
Chem. Phys. 127 , 154103 (2007). 

[33] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Phys. 
Rev. A 77 , 033613 (2008). 

[34] J. I. Kim, V. S. Melezhik, and P. Schmelcher, Phys. Rev. 
Lett. 97 , 193203 (2006). 

[35] P. Giannakeas, F. K. Diakonos, and P. Schmelcher, Phys. 
Rev. A 86, 042703 (2012). 

[36] S. I. Mistakidis, L. Cao, and P. Schmelcher, J. Phys. B: 
At. Mol. Opt. Phys. 47 , 225303 (2014). 

[37] S. I. Mistakidis, L. Cao, and P. Schmelcher, Phys. Rev. 
A 91 , 033611 (2014). 

[38] T. Gorin, T. Prosen, T. H. Seligman, and M. Znidaric, 
Phys. Rep. 435 , 33 (2006). 

[39] D.J. Tannor, Introduction to Quantum Mechanics: A 
Time-Dependent Perspective , (University Science Books, 
Sausalito, California, 2007). 


[40] M. Grifoni and P. Hanggi, Phys. Rep. 304 , 229-354 
(1998). 

[41] T. Wulf, C. Petri, B. Liebchen and P. Schmelcher, Phys. 
Rev. E 90 , 042913 (2014). 

[42] T. Wulf, B. Liebchen and P. Schmelcher, Phys. Rev. A 
91 , 043628 (2015). 

[43] L. Cao, S. Kronke, O. Vendrell, and P. Schmelcher, J. 
Chem. Phys. 139 , 134103 (2013). 

[44] S. Kronke, L. Cao, O. Vendrell, and P. Schmelcher, New 
J. Phys. 15 , 063018 (2013). 

[45] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Phys. 
Rev. Lett. 99 , 030402 (2007). 

[46] A. I. Streltsov, K. Sakmann, O. E. Alon, and L. S. Ceder¬ 
baum, Phys. Rev. A 83, 043604 (2011). 

[47] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Phys. 
Rev. A 76 , 013611 (2007). 

[48] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Phys. 
Rev. A 79 , 022503 (2009). 

[49] J. Frenkel, in Wave Mechanics 1st ed. (Clarendon Press, 
Oxford, 1934), pp. 423-428. 

[50] P. A. Dirac, (1930, July). Proc. Camb. Phil. Soc. (Vol. 
26 , No. 03, pp. 376-385). Cambridge University Press. 

[51] J. Broeckhove, L. Lathouwers, E. Kesteloot, and P. Van 
Leuven, Chem. Phys. Lett. 149 , 547 (1988). 


